#!/usr/bin/perl -w
use strict;
use Getopt::Long;

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

my $usage = <<'USAGE';

USAGE:

    create_evolver_params.pl SEQUENCE_FILES
                            [--kappa 2]
                            [--tree "(1:0.1,2:0.1)"]
                            [--omega 0.1]
                            [--cnt_output_sequences 2]
							[--replicates 8000]
                            [--length 50]
                            [--verbose]
                            [--help]
    
    Accepts a sequence file, calculates codon probabilities and
    outputs a control file for PAML evolver
    
USAGE


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory

# optional parameters
my $verbose;
my $kappa = 2;
my $omega = 0.1;
my $tree = "(1:0.1,2:0.1)";
#my $tree = "((sequence_a:0.1,sequence_b:0.1):0.1, sequence_c:0.1);";
my $length = 50;
my $replicates = 8000;
my $sequence_file;
my $help;
my $cnt_output_sequences = 2;
GetOptions(
			'verbose'					=> \$verbose,
			'sequence_file=s'			=> \$sequence_file,
			'kappa=s'					=> \$kappa,
			'omega=s'					=> \$omega,
			'tree=s'					=> \$tree,
			'length=s'					=> \$length,
			'cnt_output_sequences=s'	=> \$cnt_output_sequences,
			'replicates=s'				=> \$replicates,
			'help'						=> \$help);


die $usage if ($help);

# if tree is a floating point number
if ($tree =~
		/^
            [+-]?\ *		# first, match an optional sign *and space*
            (				# then match integers or f.p. mantissas:
				\d+\.\d+	# mantissa of the form a.b
				|\d+\.  	# mantissa of the form a.
				|\.\d+  	# mantissa of the form .b
				|\d+    	# integer of the form a
            )
            ([eE][+-]?\d+)? # finally, optionally match an exponent
		$/x)
{
	$tree *= 0.5;
	$tree = "(1:$tree,2:$tree)";
}





#
#	create array of codons in order
#
my @codon_order;
for my $pos1('T', 'C', 'A', 'G')
{
	for my $pos2('T', 'C', 'A', 'G')
	{
		for my $pos3('T', 'C', 'A', 'G')
		{
			push(@codon_order, $pos1.$pos2.$pos3);
		}
	}
}



#
#	codon frequencies
#
my %codon_freq;
@codon_freq{@codon_order} = (0) x @codon_order;


#
#	codon frequencies
#
sub accumulate_codon_frequencies($\%$)
{
	my ($seq, $codon_freq, $cnt_lines) = @_;
	die $cnt_lines unless defined $seq;
	die "$seq, $cnt_lines" unless length($seq) % 3 == 0;
	my @codons = unpack('a3'x(length($seq) / 3), $seq);
	for my $codon(@codons)
	{
		$codon_freq{$codon}++
	}
}

#
#	process sequence files
#

# initialize dots
my $dots_total = 2000000;
my $dots_expect_lines = 5;
my $dots_counter = 0;
my $dots_every = $dots_total / ($dots_expect_lines * 50) + 1;
print STDERR "\t[" if $verbose;
my $seq = '';
while (<>)
{
	# Update dots
	if (++$dots_counter % $dots_every == 0) 
	{
		my $dots_line = int($dots_counter / $dots_every);
		print STDERR "\n\t" if ($dots_line && $dots_line % 50 == 0 && $verbose); 
		print STDERR "." if $verbose;
	}
	if (/^>/)
	{
		accumulate_codon_frequencies($seq, %codon_freq, $dots_counter);
		$seq = '';
		next;
	}
	
	chomp;
	$seq .= $_;
}
# last sequence
accumulate_codon_frequencies($seq, %codon_freq, $dots_counter);
print STDERR "]\n" if $verbose;





#for (@codon_order)
#{
#	print "$_ =>$codon_freq{$_};\n";
#}

#
#	count total number of codons
#
my $cnt_total = 0;
$cnt_total += $codon_freq{$_} for (@codon_order);

#
#	print file preamble
#
print <<"CTL_FILE";
0           * 0:paml format (mc.paml); 1:paup format (mc.paup)
13147       * random number seed (odd number)
$cnt_output_sequences $length $replicates   * <# seqs>  <# nucleotide sites>  <# replicates>
-1          * <tree length; see note below; use -1 if tree has absolute branch lengths>
$tree

$omega      * omega
$kappa      * kappa
CTL_FILE


#
#	print frequencies
#
my $cnt_codons = 0;
for (@codon_order)
{
	print "\n" if ($cnt_codons % 4 == 0);
	$cnt_codons++;
	printf "%.8f ",  $codon_freq{$_} / $cnt_total;
}
print "\n// end of file.\n";

